home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / eulisp / feel0_89.lha / Feel / Modules / random.em < prev    next >
Encoding:
Text File  |  1993-07-12  |  7.6 KB  |  214 lines
  1. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  2. ;;                                                                           ;;
  3. ;;  EuLisp Module                     Copyright (C) University of Bath 1991  ;;
  4. ;;                                                                           ;;
  5. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  6.  
  7. ;;; -*- Mode: LISP; Package: random; Syntax: Common-lisp; Base: 10.;  -*- 
  8.  
  9. ;;;%Header
  10. ;;;----------------------------------------------------------------------------
  11. ;;;
  12. ;;; Pseudo-random number generator
  13. ;;;
  14. ;;; Author: Chris McConnell, ccm@cs.cmu.edu
  15. ;;;
  16. ;;; This file implements a portable pseudo-random number generator for
  17. ;;; Common LISP.  It has been converted from a C program that was
  18. ;;; converted from a FORTRAN program.  I did not pick the variable
  19. ;;; names or pretend to have figured out how it works.  The
  20. ;;; correctness of the generator can be verified by the TEST function
  21. ;;; at the end of the file.
  22. ;;;
  23. ;;; Original C header:
  24. ;;;
  25. ;;;    This is the random number generator proposed by George Marsaglia in
  26. ;;;    Florida State University Report: FSU-SCRI-87-50
  27. ;;;
  28. ;;;    This Random Number Generator is based on the algorithm in a FORTRAN
  29. ;;;    version published by George Marsaglia and Arif Zaman, Florida State
  30. ;;;    University; ref.: see original comments below.
  31. ;;; 
  32. ;;;    At the fhw (Fachhochschule Wiesbaden, W.Germany), Dept. of Computer
  33. ;;;    Science, we have written sources in further languages (C, Modula-2
  34. ;;;    Turbo-Pascal(3.0, 5.0), Basic and Ada) to get exactly the same test
  35. ;;;    results compared with the original FORTRAN version.
  36. ;;;                                                          April 1989
  37. ;;; 
  38. ;;;                                         Karl-L. Noell <NOELL@DWIFH1.BITNET>
  39. ;;;                                    and  Helmut  Weber <WEBER@DWIFH1.BITNET>
  40. ;;;
  41. ;;; Eulisp'ed: PAB
  42.  
  43. (defmodule random 
  44.   (eulisp0 loops) ()
  45.   ()
  46. ;; Extra define
  47. (defconstant mod remainder)
  48. (defun trunc (a . b)
  49.   (cond ((null b) (floor a))
  50.     (t (floor (/ a (car b))))))
  51. (defmethod floor ((x <integer>)) x)
  52. (defun minusp (a)
  53.   (< a 0))
  54.  
  55. ;;;%globals
  56. (deflocal *state* nil) ;; "the default state structure."
  57.  
  58. ;;;%%state 
  59. (defstruct state ()
  60.   ;;  "this contains random state for a state.  the names of slots are the
  61.   ;; same as the variable names in the orginal program."
  62.   ((ij initform 1802 initarg ij accessor state-ij)
  63.    (kl initform 9373 initarg kl accessor state-kl)
  64.    (u initform (make-vector 97 0) initarg u accessor state-u)
  65.    (c initform (/ 362436.0 16777216.0) initarg c accessor state-c)
  66.    (cd initform (/ 7654321.0 16777216.0) initarg cd accessor state-cd)
  67.    (cm initform (/ 16777213.0 16777216.0) accessor state-cm)
  68.    (i97 initform 96 accessor state-i97)
  69.    (j97 initform 32 accessor state-j97))
  70.   constructor make-state)
  71.  
  72. ;;;
  73. (defun print-state (state stream level)
  74.   (declare (ignore level))
  75.   (format stream "#<state ~a ~a>" (state-ij state) (state-kl state)))
  76.  
  77. ;;; %interface
  78. ;;; this is the initialization routine for the random number generator  
  79. ;;; note: the seed variables can have values between:    0 <= ij <= 31328
  80. ;;;                                                      0 <= kl <= 30081
  81. ;;; the random number sequences created by these two seeds are of sufficient
  82. ;;; length to complete an entire calculation with. for example, if sveral
  83. ;;; different groups are working on different parts of the same calculation,
  84. ;;; each group could be assigned its own ij seed. this would leave each group
  85. ;;; with 30000 choices for the second seed. that is to say, this random
  86. ;;; number generator can create 900 million different subsequences -- with
  87. ;;; each subsequence having a length of approximately 10^30.
  88. ;;; 
  89. ;; (defun seed-state (&optional (ij (mod (get-internal-real-time) 31329))
  90. ;;                 (kl (mod (get-internal-real-time) 30082))))
  91. (defun seed-state (ij kl)
  92.   ;;  "given the seed values 0 <= ij <= 31328 and 0 <= kl <= 30081,
  93.   ;; generate a state structure, set it as the default state and return it."
  94.   ;; (declare (optimize (speed 3) (safety 0)))
  95.   (let* ((state (make-state 'ij ij 'kl kl))
  96.      (vector (state-u state))
  97.      (i (+ (mod (trunc ij 177) 177) 2))
  98.      (j (+ (mod ij 177) 2))
  99.      (k (+ (mod (trunc kl 169) 178) 1))
  100.      (l (mod kl 169))
  101.      (ii 0))
  102.     ;;(declare (type fixnum i j k l)
  103.     ;; (type (simple-vector 97) vector))
  104.     (while (< ii 97)
  105.       (let ((s 0.0)
  106.         (t1 0.5)
  107.         (jj 0))
  108.     ;;(declare (type single-float s t1))
  109.     (while (< jj 24)
  110.          (let ((m (mod (* (mod (* i j) 179) k) 179)))
  111.            ;; (declare (type fixnum m))
  112.            (setq i j)
  113.            (setq j k)
  114.            (setq k m)
  115.            (setq l (mod (+ (* 53 l) 1) 169))
  116.            (if  (not (< (mod (* l m) 64) 32))
  117.            (setq s (+ s t1))
  118.          nil)
  119.            (setq t1 (* 0.5 t1)))
  120.          (setq jj (+ jj 1)))
  121.     ((setter vector-ref) vector ii s))
  122.       (setq ii (+ ii 1)))
  123.     (setq *state* state)))
  124.  
  125. ;;;
  126. (defun random-one (state)
  127.   ;; "return a random value between 0.0 and 1.0 from state."
  128.   ;; (declare (optimize (speed 3) (safety 0)))
  129.   (let* ((u (state-u state))
  130.      (uni (- (vector-ref u (state-i97 state))
  131.          (vector-ref u (state-j97 state)))))
  132.     ;;(declare (type simple-vector u)
  133.     ;;         (type single-float uni))
  134.     (when (minusp uni)
  135.       (setq uni (+ uni 1)))
  136.     ((setter vector-ref) u (state-i97 state) uni)
  137.     ((setter state-i97) state (- (state-i97 state) 1))
  138.     (when (minusp (state-i97 state))
  139.       ((setter state-i97) state 96))
  140.     ((setter state-j97) state (- (state-j97 state) 1))
  141.     (when (minusp (state-j97 state))
  142.       ((setter state-j97) state 96))
  143.     ((setter state-c) state
  144.      (- (state-c state) (state-cd state)))
  145.     (when (minusp (state-c state))
  146.       ((setter state-c) state
  147.        (+ (state-c state) (state-cm state))))
  148.     (setq uni (- uni (state-c state)))
  149.     (if (minusp uni)
  150.     (setq uni (+ uni 1))
  151.       nil)
  152.     uni))
  153.  
  154. ;;;
  155. (defun random-float (n state)
  156.   ;;"return a random float between 0 and n.
  157.   ;; use seed-state to initialize a pseudo-random sequence."
  158.   ;; (declare (optimize (speed 3) (safety 0)))
  159.   (* n (random-one state)))
  160.  
  161. ;;;
  162. (defun random-integer (n &optional (state *state*))
  163.   ;; "return a random integer between 0 and n.
  164.   ;; use seed-state to initialize a pseudo-random sequence."
  165.   (declare (optimize (speed 3) (safety 0)))
  166.   (trunc (* n (random-one state))))
  167.  
  168. ;;;
  169. (defun random-number (n &optional (state *state*))
  170.   ;; "return a random number between 0 and n with the same type as n.
  171.   ;; use seed-state to initialize a pseudo-random sequence."
  172.   ;; (declare (optimize (speed 3) (safety 0)))
  173.   (if (floatp n)
  174.       (random-float n state)
  175.       (random-integer n state)))
  176.               
  177. ;;;
  178. (defun random-range (min max &optional (state *state*))
  179.   ;; "return a number between min and max with the same type.
  180.   ;; use seed-state to initialize a pseudo-random sequence."
  181.   (declare (optimize (speed 3) (safety 0)))
  182.   (let ((number (+ min (* (- max min) (random-one state)))))
  183.     (if (floatp min)
  184.     number
  185.     (trunc number))))
  186.    
  187. ;;;%test code
  188. (defun test ()
  189.   ;; "test the random number generator.  it should print out n = n t six
  190.   ;; times if it works correctly."
  191.   (let ((state (seed-state 1802 9373))
  192.     (x 0))
  193.     (while (< x 20000)
  194.       (format t "~a: ~a~% " x (random-one state))
  195.       (setq x (+ 1 x)))
  196.     (setq x 0)
  197.     (while (< x 6)
  198.       (let ((value (round (* 4096 4096 (random-one state))))
  199.         (should-be (vector-ref #(6533892 14220222 7275067
  200.                          6172232 8354498 10633180) x)))
  201.     (format t "~&~a = ~a ~a"
  202.         value should-be (= value should-be)))
  203.       (setq x (+ 1 x)))))
  204.  
  205. ;;;%initial seed -- numbers chosen by pab.
  206. ;;(setq *state* (seed-state 23465 987458))
  207.  
  208. ;; export stuff
  209. (export *state*
  210.     seed-state random-one random-float random-integer random-range)
  211.  
  212. ;; close module
  213. )
  214.